/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package TweetRay.Backgrounds;

import Raytracing.Constants;
import processing.core.PVector;

/**
 *
 * @author Matt
 */
public class ClearSkyBackground extends Background {

    private float turbidity = 3;
    //public PVector toSun;
    /**
     * Angle in radians from the zenith to the sun vector.
     */
    public float thetaS;
    /**
     * Angle in radians to the sun vector from the south pole.
     */
    public float phiS;

    /**
     * Return the color of the sky for an angle in the sky.
     * @param dir Sky direction to sample
     * @return int color of the sky. Will be tonemapped using an exponential operator
     */
    public ClearSkyBackground() {
        thetaS = 0;
        phiS = .25f;
    }

    public void getColor(PVector dir, PVector ret) {

        if (dir.y >= 0 ){
            ret.set(0,0,0);
            return;
        }
        float theta2, theta3, T, T2, chi, zenith_Y, zenith_x, zenith_y;
        float[] perez_Y = new float[5];
        float[] perez_x = new float[5];
        float[] perez_y = new float[5];

        theta2 = thetaS * thetaS;
        theta3 = theta2 * thetaS;
        T = turbidity;
        T2 = turbidity * turbidity;

        chi = (4.f / 9.f - T / 120.f) * (Constants.PI - 2.f * thetaS);

        zenith_Y = (4.0453f * T - 4.9710f) * (float) Math.tan(chi) - .2155f * T + 2.4192f;
        zenith_Y *= 1000.f;  // conversion from kcd/m^2 to cd/m^2

        zenith_x =
                (0.00165f * theta3 - 0.00375f * theta2 + 0.00209f * thetaS + 0.0f) * T2 +
                (-0.02903f * theta3 + 0.06377f * theta2 - 0.03202f * thetaS + 0.00394f) * T +
                (0.11693f * theta3 - 0.21196f * theta2 + 0.06052f * thetaS + 0.25886f);

        zenith_y =
                (0.00275f * theta3 - 0.00610f * theta2 + 0.00317f * thetaS + 0.0f) * T2 +
                (-0.04214f * theta3 + 0.08970f * theta2 - 0.04153f * thetaS + 0.00516f) * T +
                (0.15346f * theta3 - 0.26756f * theta2 + 0.06670f * thetaS + 0.26688f);

        perez_Y[0] = 0.17872f * T - 1.46303f;
        perez_Y[1] = -0.35540f * T + 0.42749f;
        perez_Y[2] = -0.02266f * T + 5.32505f;
        perez_Y[3] = 0.12064f * T - 2.57705f;
        perez_Y[4] = -0.06696f * T + 0.37027f;

        perez_x[0] = -0.01925f * T - 0.25922f;
        perez_x[1] = -0.06651f * T + 0.00081f;
        perez_x[2] = -0.00041f * T + 0.21247f;
        perez_x[3] = -0.06409f * T - 0.89887f;
        perez_x[4] = -0.00325f * T + 0.04517f;

        perez_y[0] = -0.01669f * T - 0.26078f;
        perez_y[1] = -0.09495f * T + 0.00921f;
        perez_y[2] = -0.00792f * T + 0.21023f;
        perez_y[3] = -0.04405f * T - 1.65369f;
        perez_y[4] = -0.01092f * T + 0.05291f;


        float theta = (float) Math.acos(-dir.y);
        float phi = (float) Math.atan2(dir.z, dir.x);



        float gamma = AngleBetween(theta, phi, thetaS, phiS);

        // Compute xyY values
        float x = perezFunction(perez_x, gamma, theta, thetaS, zenith_x);
        float y = perezFunction(perez_y, gamma, theta, thetaS, zenith_y);
        float Y = perezFunction(perez_Y, gamma, theta, thetaS, zenith_Y);

        float X = (x / y) * Y;
        float Z = ((1.f - x - y) / y) * Y;


        float[] rgb = new float[3];
        rgb[0] = 1.f - (float) Math.exp(-(1.f / 15000.f) * (3.240479f * X - 1.537150f * Y - 0.498535f * Z));
        rgb[1] = 1.f - (float) Math.exp(-(1.f / 15000.f) * (-0.969256f * X + 1.875992f * Y + 0.041556f * Z));
        rgb[2] = 1.f - (float) Math.exp(-(1.f / 15000.f) * (0.055648f * X - 0.204043f * Y + 1.057311f * Z));

        if (rgb[0] < 0) rgb[0] = 0;
        if (rgb[1] < 0) rgb[1] = 0;
        if (rgb[2] < 0) rgb[2] = 0;
        ret.set(rgb[0], rgb[1], rgb[2]);



    }

    private float AngleBetween(float thetav, float phiv, float theta, float phi) {
        float cospsi = (float) (Math.sin(thetav) * Math.sin(theta) * Math.cos(phi - phiv) + Math.cos(thetav) * Math.cos(theta));
        if (cospsi > 1) {
            return 0;
        }
        if (cospsi < -1) {
            return Constants.PI;
        }
        return (float) Math.acos(cospsi);
    }

    private float perezFunction(float[] coeffs, float gamma, float theta, float thetaS, float lvz) {
        float A = coeffs[0], B = coeffs[1], C = coeffs[2], D = coeffs[3], E = coeffs[4];

        float den = (float) ((1.f + A * Math.exp(B)) * (1.f + C * Math.exp(D * thetaS) + E * Math.cos(thetaS) * Math.cos(thetaS)));
        float num = (float) ((1.f + A * Math.exp(B / Math.cos(theta))) * (1.f + C * Math.exp(D * gamma) + E * Math.cos(gamma) * Math.cos(gamma)));
        return (lvz * num / den);

    }



    private float radians(float x) {
        return x / 180.f * Constants.PI;

    }
}
